/*
 * com.javable.cese -- Cell Electrophysiology Simulation Environment (CESE)
 *
 * Copyright (C) 1999 - 2005 Sergey Missan
 *
 * Dalhousie Univ., Halifax, NS, Canada
 *
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public
 * License along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
 *
 */

package models.com.javable.cese.models;

import java.io.Serializable;

import com.javable.cese.templates.AbstractAPModel;

/**
 * The Updated Luo-Rudy II Model of Mammalian Ventricular Cardiac Action
 * Potentials, 2000. Detailed list of equations and model description are
 * provided in
 * <p>
 * Circ Res 1991;68:1501-1526 <br>
 * Circ Res 1994;74:1071-1096 <br>
 * Circ Res 1994;74:1097-1113 <br>
 * Circ Res 1995;77:140-152 <br>
 * Biophys J 1995;68:949-964 <br>
 * Cardiovasc Res 1997;35:256-272 <br>
 * Circulation 1999;99:2466-2474 <br>
 * Cardiovas Res 1999;42:530-542 <br>
 * Nature 1999;400:566-569 <br>
 * Circulation 2000;101:1192-1198 <br>
 * Biophy J 2000;78:2392-2404.
 */
public class LuoRudyModel_II_2000 extends AbstractAPModel implements
        Serializable {

    /* The beginning here lists all the variables and defines them */

    /* Cell Geometry */
    private double l = 0.01; // Length of the cell (cm)

    public synchronized void setl(double sl) {
        l = sl;
        getCellGeometry();
    }

    public double getl() {
        return l;
    }

    private double a = 0.0011; // Radius of the cell (cm)

    public synchronized void seta(double sa) {
        a = sa;
        getCellGeometry();
    }

    public double geta() {
        return a;
    }

    private double vcell; // Cell volume (uL)

    public double getvcell() {
        return vcell;
    }

    private double ageo; // Geometric membrane area (cm^2)

    public double getageo() {
        return ageo;
    }

    private double acap; // Capacitive membrane area (cm^2)

    public double getacap() {
        return acap;
    }

    private double vmyo; // Myoplasm volume (uL)

    public double getvmyo() {
        return vmyo;
    }

    private double vmito; // Mitochondria volume (uL)

    public double getvmito() {
        return vmito;
    }

    private double vsr; // SR volume (uL)

    public double getvsr() {
        return vsr;
    }

    private double vnsr; // NSR volume (uL)

    public double getvnsr() {
        return vnsr;
    }

    private double vjsr; // JSR volume (uL)

    public double getvjsr() {
        return vjsr;
    }

    private double vcleft; // Cleft volume (uL)

    public double getvcleft() {
        return vcleft;
    }

    /* Terms for Solution of Conductance and Reversal Potential */
    public static final double R = 8314; // Universal Gas Constant (J/kmol*K)
    public static final double frdy = 96485; // Faraday's Constant (C/mol)

    private double temp = 310; // Temperature (K)

    public synchronized void settemp(double stemp) {
        temp = stemp;
    }

    public double gettemp() {
        return temp;
    }

    /* Ion Valences */
    public static final double zna = 1; // Na valence
    public static final double zk = 1; // K valence
    public static final double zca = 2; // Ca valence

    /* Ion Concentrations */
    private double nai = 10; // Intracellular Na Concentration (mM)

    public synchronized void setnai(double snai) {
        nai = snai;
    }

    public double getnai() {
        return nai;
    }

    private double nao = 140; // Extracellular Na Concentration (mM)

    public synchronized void setnao(double snao) {
        nao = snao;
    }

    public double getnao() {
        return nao;
    }

    private double ki = 145; // Intracellular K Concentration (mM)

    public synchronized void setki(double ski) {
        ki = ski;
    }

    public double getki() {
        return ki;
    }

    private double ko = 5.4; // Extracellular K Concentration (mM)

    public synchronized void setko(double sko) {
        ko = sko;
    }

    public double getko() {
        return ko;
    }

    private double cai = 0.00012; // Intracellular Ca Concentration (mM)

    public synchronized void setcai(double scai) {
        cai = scai;
    }

    public double getcai() {
        return cai;
    }

    private double cao = 1.8; // Extracellular Ca Concentration (mM)

    public synchronized void setcao(double scao) {
        cao = scao;
    }

    public double getcao() {
        return cao;
    }

    private double cmdn = 0.0024; // Calmodulin Buffered Ca Concentration (mM)

    public synchronized void setcmdn(double scmdn) {
        cmdn = scmdn;
    }

    public double getcmdn() {
        return cmdn;
    }

    private double trpn = 0.0133; // Troponin Buffered Ca Concentration (mM)

    public synchronized void settrpn(double strpn) {
        trpn = strpn;
    }

    public double gettrpn() {
        return trpn;
    }

    private double nsr = 1.8; // NSR Ca Concentration (mM)

    public synchronized void setnsr(double snsr) {
        nsr = snsr;
    }

    public double getnsr() {
        return nsr;
    }

    private double jsr = 1.8; // JSR Ca Concentration (mM)

    public synchronized void setjsr(double sjsr) {
        jsr = sjsr;
    }

    public double getjsr() {
        return jsr;
    }

    private double csqn = 6.9; // Calsequestrin Buffered Ca Concentration (mM)

    public synchronized void setcsqn(double scsqn) {
        csqn = scsqn;
    }

    public double getcsqn() {
        return csqn;
    }

    /* Myoplasmic Na Ion Concentration Changes */
    private double naiont; // Total Na Ion Flow (mM/ms)

    public double getnaiont() {
        return naiont;
    }

    private double dnai; // Change in Intracellular Na Concentration (mM)

    public double getdnai() {
        return dnai;
    }

    /* Myoplasmic K Ion Concentration Changes */
    private double kiont; // Total K Ion Flow (mM/ms)

    public double getkiont() {
        return kiont;
    }

    private double dki; // Change in Intracellular K Concentration (mM)

    public double getdki() {
        return dki;
    }

    /* NSR Ca Ion Concentration Changes */
    private double dnsr; // Change in [Ca] in the NSR (mM)

    public double getdnsr() {
        return dnsr;
    }

    private double iup; // Ca uptake from myo. to NSR (mM/ms)

    public double getiup() {
        return iup;
    }

    private double ileak; // Ca leakage from NSR to myo. (mM/ms)

    public double getileak() {
        return ileak;
    }

    private double kleak; // Rate constant of Ca leakage from NSR (ms^-1)

    public double getkleak() {
        return kleak;
    }

    private double kmup = 0.00092; // Half-saturation concentration of iup (mM)

    public synchronized void setkmup(double skmup) {
        kmup = skmup;
    }

    public double getkmup() {
        return kmup;
    }

    private double iupbar = 0.005; // Max. current through iup channel (mM/ms)

    public synchronized void setiupbar(double siupbar) {
        iupbar = siupbar;
    }

    public double getiupbar() {
        return iupbar;
    }

    private double nsrbar = 15; // Max. [Ca] in NSR (mM)

    public synchronized void setnsrbar(double snsrbar) {
        nsrbar = snsrbar;
    }

    public double getnsrbar() {
        return nsrbar;
    }

    /* JSR Ca Ion Concentration Changes */
    private double djsr; // Change in [Ca] in the JSR (mM)

    public double getdjsr() {
        return djsr;
    }

    private double tauon = 2;

    // Time constant of activation of Ca release from JSR (ms)
    public synchronized void settauon(double stauon) {
        tauon = stauon;
    }

    public double gettauon() {
        return tauon;
    }

    private double tauoff = 2;

    // Time constant of deactivation of Ca release from JSR (ms)
    public synchronized void settauoff(double stauoff) {
        tauoff = stauoff;
    }

    public double gettauoff() {
        return tauoff;
    }

    private double caitzero;

    // Total myoplasmic [Ca] at time of dv/dt max. (mM)
    public double getcaitzero() {
        return caitzero;
    }

    private double caiatwo;

    // Total myoplasmic [Ca] 2 ms. after dv/dt max. (mM)
    public double getcaiatwo() {
        return caiatwo;
    }

    private double dcaitwo;

    // Cumulative Ca entry into the cell 2 ms. after dv/dt max. (mM)
    public double getdcaitwo() {
        return dcaitwo;
    }

    private double tcount = 3;
    // Counter for 2 ms. after onset of dv/dt max. (ms)

    private double dcaith = 0.000242;

    // Threshold for external triggering of Ca release from JSR (mM)
    public synchronized void setdcaith(double sdcaith) {
        dcaith = sdcaith;
    }

    public double getdcaith() {
        return dcaith;
    }

    private double grelbarcicr = 0;

    // Max. rate constant of Ca release from JSR due to CICR (ms^-1)
    public synchronized void setgrelbarcicr(double sgrelbarcicr) {
        grelbarcicr = sgrelbarcicr;
    }

    public double getgrelbarcicr() {
        return grelbarcicr;
    }

    private double grelcicr;

    // Rate constant of Ca release from JSR due to CICR (ms^-1)
    public double getgrelcicr() {
        return grelcicr;
    }

    private double tcicr; // t=0 at time of CICR (ms)

    public double gettcicr() {
        return tcicr;
    }

    private double kmrel = 0.0008;

    // Half-saturation constant of irel for CICR (mM)
    public synchronized void setkmrel(double skmrel) {
        kmrel = skmrel;
    }

    public double getkmrel() {
        return kmrel;
    }

    private double irelcicr; // Ca release from JSR to myo. due to CICR (mM/ms)

    public double getirelcicr() {
        return irelcicr;
    }

    private double csqnth = 8.75;

    // Threshold for release of Ca from Calsequestrin (CSQN) due to JSR
    // overload (mM)
    public synchronized void setcsqnth(double scsqnth) {
        csqnth = scsqnth;
    }

    public double getcsqnth() {
        return csqnth;
    }

    private double gmaxrel = 60;

    // Max. rate constant of Ca release from JSR due to overload (ms^-1)
    public synchronized void setgmaxrel(double sgmaxrel) {
        gmaxrel = sgmaxrel;
    }

    public double getgmaxrel() {
        return gmaxrel;
    }

    private double grelbarjsrol = 0;

    // Rate constant of Ca release from JSR due to overload (ms^-1)
    public synchronized void setgrelbarjsrol(double sgrelbarjsrol) {
        grelbarjsrol = sgrelbarjsrol;
    }

    public double getgrelbarjsrol() {
        return grelbarjsrol;
    }

    private double greljsrol;

    // Rate constant of Ca release from JSR due to JSR overload (ms^-1)
    public double getgreljsrol() {
        return greljsrol;
    }

    private double tjsrol; // t=0 at time of JSR overload (ms)

    public double gettjsrol() {
        return tjsrol;
    }

    private double ireljsrol;

    // Ca release from JSR to myo. due to JSR overload (mM/ms)
    public double getireljsrol() {
        return ireljsrol;
    }

    private double csqnbar = 10; // Max. [Ca] buffered in CSQN (mM)

    public synchronized void setcsqnbar(double scsqnbar) {
        csqnbar = scsqnbar;
    }

    public double getcsqnbar() {
        return csqnbar;
    }

    private double kmcsqn = 0.8;

    // Equalibrium constant of buffering for CSQN (mM)
    public synchronized void setkmcsqn(double skmcsqn) {
        kmcsqn = skmcsqn;
    }

    public double getkmcsqn() {
        return kmcsqn;
    }

    private double bjsr;

    // b Variable for analytical computation of [Ca] in JSR (mM)
    public double getbjsr() {
        return bjsr;
    }

    private double cjsr;

    // c Variable for analytical computation of [Ca] in JSR (mM)
    public double getcjsr() {
        return cjsr;
    }

    /* Translocation of Ca Ions from NSR to JSR */
    private double itr;

    // Translocation current of Ca ions from NSR to JSR (mM/ms)
    public double getitr() {
        return itr;
    }

    private double tautr = 180;

    // Time constant of Ca transfer from NSR to JSR (ms)
    public synchronized void settautr(double stautr) {
        tautr = stautr;
    }

    public double gettautr() {
        return tautr;
    }

    /* Myoplasmic Ca Ion Concentration Changes */
    private double caiont; // Total Ca Ion Flow (mM/ms)

    public double getcaiont() {
        return caiont;
    }

    private double dcai; // Change in myoplasmic Ca concentration (mM)

    public double getdcai() {
        return dcai;
    }

    private double catotal; // Total myoplasmic Ca concentration (mM)

    public double getcatotal() {
        return catotal;
    }

    private double bmyo;

    // b Variable for analytical computation of [Ca] in myoplasm (mM)
    public double getbmyo() {
        return bmyo;
    }

    private double cmyo;

    // c Variable for analytical computation of [Ca] in myoplasm (mM)
    public double getcmyo() {
        return cmyo;
    }

    private double dmyo;

    // d Variable for analytical computation of [Ca] in myoplasm (mM)
    public double getdmyo() {
        return dmyo;
    }

    private double cmdnbar = 0.050;

    // Max. [Ca] buffered in Calmodulin (CMDN) (mM)
    public synchronized void setcmdnbar(double scmdnbar) {
        cmdnbar = scmdnbar;
    }

    public double getcmdnbar() {
        return cmdnbar;
    }

    private double trpnbar = 0.070;

    // Max. [Ca] buffered in Troponin (TRPN) (mM)
    public synchronized void settrpnbar(double strpnbar) {
        trpnbar = strpnbar;
    }

    public double gettrpnbar() {
        return trpnbar;
    }

    private double kmcmdn = 0.00238;

    // Equalibrium constant of buffering for CMDN (mM)
    public synchronized void setkmcmdn(double skmcmdn) {
        kmcmdn = skmcmdn;
    }

    public double getkmcmdn() {
        return kmcmdn;
    }

    private double kmtrpn = 0.0005;

    // Equalibrium constant of buffering for TRPN (mM)
    public synchronized void setkmtrpn(double skmtrpn) {
        kmtrpn = skmtrpn;
    }

    public double getkmtrpn() {
        return kmtrpn;
    }

    private double gpig;
    // Tribute to all the guinea pigs killed for the advancement of knowledge

    /* Fast Sodium Current (time dependent) */
    private double ina; // Fast Na Current (uA/uF)

    public synchronized void setina(double sina) {
        ina = sina;
    }

    public double getina() {
        return ina;
    }

    private double gna; // Max. Conductance of the Na Channel (mS/uF)

    public synchronized void setgna(double sgna) {
        gna = sgna;
    }

    public double getgna() {
        return gna;
    }

    private double ena; // Reversal Potential of Na (mV)

    public synchronized void setena(double sena) {
        ena = sena;
    }

    public double getena() {
        return ena;
    }

    private double ah; // Na alpha-h rate constant (ms^-1)

    public double getah() {
        return ah;
    }

    private double bh; // Na beta-h rate constant (ms^-1)

    public double getbh() {
        return bh;
    }

    private double aj; // Na alpha-j rate constant (ms^-1)

    public double getaj() {
        return aj;
    }

    private double bj; // Na beta-j rate constant (ms^-1)

    public double getbj() {
        return bj;
    }

    private double am; // Na alpha-m rate constant (ms^-1)

    public double getam() {
        return am;
    }

    private double bm; // Na beta-m rate constant (ms^-1)

    public double gebm() {
        return bm;
    }

    private double h; // Na inactivation h

    public double geth() {
        return h;
    }

    private double j; // Na inactivation j

    public double getj() {
        return j;
    }

    private double m; // Na activation m

    public double getm() {
        return m;
    }

    private double scgna = 1;

    // Scaling factor for Max. Conductance of the Na Channel
    public synchronized void setscgna(double sscgna) {
        scgna = sscgna;
    }

    public double getscgna() {
        return scgna;
    }

    /* Current through L-type Ca Channel */
    private double ilca; // Ca current through L-type Ca channel (uA/uF)

    public synchronized void setilca(double silca) {
        ilca = silca;
    }

    public double getilca() {
        return ilca;
    }

    private double ilcana; // Na current through L-type Ca channel (uA/uF)

    public synchronized void setilcana(double silcana) {
        ilcana = silcana;
    }

    public double getilcana() {
        return ilcana;
    }

    private double ilcak; // K current through L-type Ca channel (uA/uF)

    public synchronized void setilcak(double silcak) {
        ilcak = silcak;
    }

    public double getilcak() {
        return ilcak;
    }

    private double ilcatot;

    // Total current through the L-type Ca channel (uA/uF)
    public double getilcatot() {
        return ilcatot;
    }

    private double ibarca; // Max. Ca current through Ca channel (uA/uF)

    public double getibarca() {
        return ibarca;
    }

    private double ibarna; // Max. Na current through Ca channel (uA/uF)

    public double getibarna() {
        return ibarna;
    }

    private double ibark; // Max. K current through Ca channel (uA/uF)

    public double getibark() {
        return ibark;
    }

    private double ad; // Ca alpha-d rate constant (ms^-1)

    public double getad() {
        return ad;
    }

    private double bd; // Ca beta-d rate constant (ms^-1)

    public double getbd() {
        return bd;
    }

    private double d; // Voltage dependent activation gate d

    public double getd() {
        return d;
    }

    private double dss; // Steady-state value of activation gate d

    public double getdss() {
        return dss;
    }

    private double taud; // Time constant of gate d (ms^-1)

    public double gettaud() {
        return taud;
    }

    private double af; // Ca alpha-f rate constant (ms^-1)

    public double getaf() {
        return af;
    }

    private double bf; // Ca beta-f rate constant (ms^-1)

    public double getbf() {
        return bf;
    }

    private double f; // Voltage dependent inactivation gate f

    public double getf() {
        return f;
    }

    private double fss; // Steady-state value of inactivation gate f

    public double getfss() {
        return fss;
    }

    private double tauf; // Time constant of gate f (ms^-1)

    public double gettauf() {
        return tauf;
    }

    private double fca; // Ca dependent inactivation gate

    public double getfca() {
        return fca;
    }

    private double kmca = 0.0006;

    // Half-saturation concentration of Ca channel (mM)
    public synchronized void setkmca(double skmca) {
        kmca = skmca;
    }

    public double getkmca() {
        return kmca;
    }

    private double pca = 0.00054; // Permiability of membrane to Ca (cm/s)

    public synchronized void setpca(double spca) {
        pca = spca;
    }

    public double getpca() {
        return pca;
    }

    private double gacai = 1; // Activity coefficient of Ca in

    public synchronized void setgacai(double sgacai) {
        gacai = sgacai;
    }

    public double getgacai() {
        return gacai;
    }

    private double gacao = 0.341; // Activity coefficient of Ca out

    public synchronized void setgacao(double sgacao) {
        gacao = sgacao;
    }

    public double getgacao() {
        return gacao;
    }

    private double pna = 0.000000675; // Permiability of membrane to Na (cm/s)

    public synchronized void setpna(double spna) {
        pna = spna;
    }

    public double getpna() {
        return pna;
    }

    private double ganai = 0.75; // Activity coefficient of Na in

    public synchronized void setganai(double sganai) {
        ganai = sganai;
    }

    public double getganai() {
        return ganai;
    }

    private double ganao = 0.75; // Activity coefficient of Na out

    public synchronized void setganao(double sganao) {
        ganao = sganao;
    }

    public double getganao() {
        return ganao;
    }

    private double pk = 0.000000193; // Permiability of membrane to K (cm/s)

    public synchronized void setpk(double spk) {
        pk = spk;
    }

    public double getpk() {
        return pk;
    }

    private double gaki = 0.75; // Activity coefficient of K in

    public synchronized void setgaki(double sgaki) {
        gaki = sgaki;
    }

    public double getgaki() {
        return gaki;
    }

    private double gako = 0.75; // Activity coefficient of K out

    public synchronized void setgako(double sgako) {
        gako = sgako;
    }

    public double getgako() {
        return gako;
    }

    private double scgcal = 1;

    // Scaling factor for Max. Conductance of the L-type Ca channel
    public synchronized void setscgcal(double sscgcal) {
        scgcal = sscgcal;
    }

    public double getscgcal() {
        return scgcal;
    }

    /* Current through T-type Ca Channel */
    private double icat; // Ca current through T-type Ca channel (uA/uF)

    public synchronized void seticat(double sicat) {
        icat = sicat;
    }

    public double geticat() {
        return icat;
    }

    private double gcat; // Max. Conductance of the T-type Ca channel (mS/uF)

    public synchronized void setgcat(double sgcat) {
        gcat = sgcat;
    }

    public double getgcat() {
        return gcat;
    }

    private double eca; // Reversal Potential of the T-type Ca channel (mV)

    public synchronized void seteca(double seca) {
        eca = seca;
    }

    public double geteca() {
        return eca;
    }

    private double ab; // Ca alpha-b rate constant (ms^-1)

    public double getab() {
        return ab;
    }

    private double bb; // Ca beta-b rate constant (ms^-1)

    public double getbb() {
        return bb;
    }

    private double b; // Voltage dependent activation gate b

    public double getb() {
        return b;
    }

    private double bss; // Steady-state value of activation gate b

    public double getbss() {
        return bss;
    }

    private double taub; // Time constant of gate b (ms^-1)

    public double gettaub() {
        return taub;
    }

    private double ag; // Ca alpha-g rate constant (ms^-1)

    public double getag() {
        return ag;
    }

    private double bg; // Ca beta-g rate constant (ms^-1)

    public double getbg() {
        return bg;
    }

    private double g; // Voltage dependent inactivation gate g

    public double getg() {
        return g;
    }

    private double gss; // Steady-state value of inactivation gate g

    public double getgss() {
        return gss;
    }

    private double taug; // Time constant of gate g (ms^-1)

    public double gettaug() {
        return taug;
    }

    private double scgcat = 1;

    // Scaling factor for Max. Conductance of the T-type Ca channel
    public synchronized void setscgcat(double sscgcat) {
        scgcat = sscgcat;
    }

    public double getscgcat() {
        return scgcat;
    }

    /* Rapidly Activating Potassium Current */
    private double ikr; // Rapidly Activating K Current (uA/uF)

    public synchronized void setikr(double sikr) {
        ikr = sikr;
    }

    public double getikr() {
        return ikr;
    }

    private double gkr;

    // Channel Conductance of Rapidly Activating K Current (mS/uF)
    public synchronized void setgkr(double sgkr) {
        gkr = sgkr;
    }

    public double getgkr() {
        return gkr;
    }

    private double ekr;

    // Reversal Potential of Rapidly Activating K Current (mV)
    public synchronized void setekr(double sekr) {
        ekr = sekr;
    }

    public double getekr() {
        return ekr;
    }

    private double axr; // K alpha-xr rate constant (ms^-1)

    public double getaxr() {
        return axr;
    }

    private double bxr; // K beta-xr rate constant (ms^-1)

    public double getbxr() {
        return bxr;
    }

    private double xr; // Rapidly Activating K time-dependent activation

    public double getxr() {
        return xr;
    }

    private double xrss; // Steady-state value of inactivation gate xr

    public double getxrss() {
        return xrss;
    }

    private double tauxr; // Time constant of gate xr (ms^-1)

    public double gettauxr() {
        return tauxr;
    }

    private double r; // K time-independent inactivation r

    public double getr() {
        return r;
    }

    private double scgkr = 1;

    // Scale factor for Channel Conductance of Rapidly Activating K Current
    public synchronized void setscgkr(double sscgkr) {
        scgkr = sscgkr;
    }

    public double getscgkr() {
        return scgkr;
    }

    /* Slowly Activating Potassium Current */
    private double iks; // Slowly Activating K Current (uA/uF)

    public synchronized void setiks(double siks) {
        iks = siks;
    }

    public double getiks() {
        return iks;
    }

    private double gks;

    // Channel Conductance of Slowly Activating K Current (mS/uF)
    public synchronized void setgks(double sgks) {
        gks = sgks;
    }

    public double getgks() {
        return gks;
    }

    private double eks;

    // Reversal Potential of Slowly Activating K Current (mV)
    public synchronized void seteks(double seks) {
        eks = seks;
    }

    public double geteks() {
        return eks;
    }

    private double axs1; // K alpha-xs1 rate constant (ms^-1)

    public double getaxs1() {
        return axs1;
    }

    private double axs2; // K alpha-xs2 rate constant (ms^-1)

    public double getaxs2() {
        return axs2;
    }

    private double bxs1; // K beta-xs1 rate constant (ms^-1)

    public double getbxs1() {
        return bxs1;
    }

    private double bxs2; // K beta-xs2 rate constant (ms^-1)

    public double getbxs2() {
        return bxs2;
    }

    private double xs1; // Slowly Activating K time-dependent activation 1

    public double getxs1() {
        return xs1;
    }

    private double xs2; // Slowly Activating K time-dependent activation 2

    public double getxs2() {
        return xs2;
    }

    private double xsss; // Steady-state value of inactivation gate xs

    public double getxsss() {
        return xsss;
    }

    private double tauxs1; // Time constant 1 of gate xs (ms^-1)

    public double gettauxs1() {
        return tauxs1;
    }

    private double tauxs2; // Time constant 2 of gate xs (ms^-1)

    public double gettauxs2() {
        return tauxs2;
    }

    private double prnak = 0.01833; // Na/K Permiability Ratio

    public synchronized void setprnak(double sprnak) {
        prnak = sprnak;
    }

    public double getprnak() {
        return prnak;
    }

    private double scgks = 1;

    // Scaling factor for Channel Conductance of Slowly Activating K Current
    public synchronized void setscgks(double sscgks) {
        scgks = sscgks;
    }

    public double getscgks() {
        return scgks;
    }

    /* Potassium Current (time-independent) */
    private double iki; // Time-independent K current (uA/uF)

    public synchronized void setiki(double siki) {
        iki = siki;
    }

    public double getiki() {
        return iki;
    }

    private double gki;

    // Channel Conductance of Time Independent K Current (mS/uF)
    public synchronized void setgki(double sgki) {
        gki = sgki;
    }

    public double getgki() {
        return gki;
    }

    private double eki;

    // Reversal Potential of Time Independent K Current (mV)
    public synchronized void seteki(double seki) {
        eki = seki;
    }

    public double geteki() {
        return eki;
    }

    private double aki; // K alpha-ki rate constant (ms^-1)

    public double getaki() {
        return aki;
    }

    private double bki; // K beta-ki rate constant (ms^-1)

    public double getbki() {
        return bki;
    }

    private double kin; // K inactivation

    public double getkin() {
        return kin;
    }

    private double scgki = 1;

    // Scaling factor for Channel Conductance of Time Independent K Current
    public synchronized void setscgki(double sscgki) {
        scgki = sscgki;
    }

    public double getscgki() {
        return scgki;
    }

    /* Plateau Potassium Current */
    private double ikp; // Plateau K current (uA/uF)

    public synchronized void setikp(double sikp) {
        ikp = sikp;
    }

    public double getikp() {
        return ikp;
    }

    private double gkp; // Channel Conductance of Plateau K Current (mS/uF)

    public synchronized void setgkp(double sgkp) {
        gkp = sgkp;
    }

    public double getgkp() {
        return gkp;
    }

    private double ekp; // Reversal Potential of Plateau K Current (mV)

    public synchronized void setekp(double sekp) {
        ekp = sekp;
    }

    public double getekp() {
        return ekp;
    }

    private double kp; // K plateau factor

    public double getkp() {
        return kp;
    }

    private double scgkp = 1;

    // Scaling factor for Channel Conductance of Plateau K Current
    public synchronized void setscgkp(double sscgkp) {
        scgkp = sscgkp;
    }

    public double getscgkp() {
        return scgkp;
    }

    /* Sodium-Calcium Exchanger */
    private double inaca; // NaCa exchanger current (uA/uF)

    public synchronized void setinaca(double sinaca) {
        inaca = sinaca;
    }

    public double getinaca() {
        return inaca;
    }

    private double knaca = 2000; // Scaling factor for inaca (uA/uF)

    public synchronized void setknaca(double sknaca) {
        knaca = sknaca;
    }

    public double getknaca() {
        return knaca;
    }

    private double kmnae = 87.5;

    // Half-saturation concentration of NaCa exhanger Na (mM)
    public synchronized void setkmnae(double skmnae) {
        kmnae = skmnae;
    }

    public double getkmnae() {
        return kmnae;
    }

    private double kmcae = 1.38;

    // Half-saturation concentration of NaCa exhanger Ca (mM)
    public synchronized void setkmcae(double skmcae) {
        kmcae = skmcae;
    }

    public double getkmcae() {
        return kmcae;
    }

    private double ksat = 0.1;

    // Saturation factor of INaCa at very negative potentials
    public synchronized void setksat(double sksat) {
        ksat = sksat;
    }

    public double getksat() {
        return ksat;
    }

    private double nu = 0.35;

    // Position of energy barrier controlling voltage dependance of inaca
    public synchronized void setnu(double snu) {
        nu = snu;
    }

    public double getnu() {
        return nu;
    }

    /* Sodium-Potassium Pump */
    private double inak; // NaK pump current (uA/uF)

    public synchronized void setinak(double sinak) {
        inak = sinak;
    }

    public double getinak() {
        return inak;
    }

    private double fnak; // Voltage-dependence parameter of inak

    public double getfnak() {
        return fnak;
    }

    private double sigma; // [Na]o dependence factor of fnak

    public double getsigma() {
        return sigma;
    }

    private double ibarnak = 1.5; // Max. current through Na-K pump (uA/uF)

    public synchronized void setibarnak(double sibarnak) {
        ibarnak = sibarnak;
    }

    public double getibarnak() {
        return ibarnak;
    }

    private double kmnai = 10;

    // Half-saturation concentration of NaK pump Na (mM)
    public synchronized void setkmnai(double skmnai) {
        kmnai = skmnai;
    }

    public double getkmnai() {
        return kmnai;
    }

    private double kmko = 1.5;

    // Half-saturation concentration of NaK pump K (mM)
    public synchronized void setkmko(double skmko) {
        kmko = skmko;
    }

    public double getkmko() {
        return kmko;
    }

    /* Nonspecific Ca-activated Current */
    private double insna; // Non-specific Na current (uA/uF)

    public synchronized void setinsna(double sinsna) {
        insna = sinsna;
    }

    public double getinsna() {
        return insna;
    }

    private double insk; // Non-specific K current (uA/uF)

    public synchronized void setinsk(double sinsk) {
        insk = sinsk;
    }

    public double getinsk() {
        return insk;
    }

    private double ibarnsna; // Max. Na current through NSCa channel (uA/uF)

    public double getibarnsna() {
        return ibarnsna;
    }

    private double ibarnsk; // Max. K current through NSCa channel (uA/uF)

    public double getibarnsk() {
        return ibarnsk;
    }

    private double pnsca = 0.000000175;

    // Permiability of membrane to Ca (cm/s)
    public synchronized void setpnsca(double spnsca) {
        pnsca = spnsca;
    }

    public double getpnsca() {
        return pnsca;
    }

    private double kmnsca = 0.0012;

    // Half-saturation concentration of NSCa channel (mM)
    public synchronized void setkmnsca(double skmnsca) {
        kmnsca = skmnsca;
    }

    public double getkmnsca() {
        return kmnsca;
    }

    /* ATP-sensitive Potassium Current */
    protected double mgi = 1.0; // Intracellular Mg Concentration (mM)

    public synchronized void setmgi(double smgi) {
        mgi = smgi;
    }

    public double getmgi() {
        return mgi;
    }

    protected double atpi = 6800; // Intracellular ATP Concentration (uM)

    public synchronized void setatpi(double satpi) {
        atpi = satpi;
    }

    public double getatpi() {
        return atpi;
    }

    protected double adpi = 15; // Intracellular ADP Concentration (uM)

    public synchronized void setadpi(double sadpi) {
        adpi = sadpi;
    }

    public double getadpi() {
        return adpi;
    }

    private double ikatp; // ATP-sensitive Potassium Current (uA/uF)

    public synchronized void setikatp(double sikatp) {
        ikatp = sikatp;
    }

    public double getikatp() {
        return ikatp;
    }

    private double fatp; // Fraction of open KATP channels

    public double getfatp() {
        return fatp;
    }

    private double gbarkatp; // IKATP Conductance

    public double getgbarkatp() {
        return gbarkatp;
    }

    private double scgbarkatp = 0; // Scaling factor for IKATP Conductance

    public synchronized void setscgbarkatp(double sscgbarkatp) {
        scgbarkatp = sscgbarkatp;
    }

    public double getscgbarkatp() {
        return scgbarkatp;
    }

    /* Sarcolemmal Ca Pump */
    private double ipca; // Sarcolemmal Ca pump current (uA/uF)

    public synchronized void setipca(double sipca) {
        ipca = sipca;
    }

    public double getipca() {
        return ipca;
    }

    private double ibarpca = 1.15;

    // Max. Ca current through sarcolemmal Ca pump (uA/uF)
    public synchronized void setibarpca(double sibarpca) {
        ibarpca = sibarpca;
    }

    public double getibarpca() {
        return ibarpca;
    }

    private double kmpca = 0.0005;

    // Half-saturation concentration of sarcolemmal Ca pump (mM)
    public synchronized void setkmpca(double skmpca) {
        kmpca = skmpca;
    }

    public double getkmpca() {
        return kmpca;
    }

    /* Ca Background Current */
    private double icab; // Ca background current (uA/uF)

    public synchronized void seticab(double sicab) {
        icab = sicab;
    }

    public double geticab() {
        return icab;
    }

    private double gcab; // Max. conductance of Ca background (mS/uF)

    public double getgcab() {
        return gcab;
    }

    private double ecan; // Nernst potential for Ca (mV)

    public synchronized void setecan(double secan) {
        ecan = secan;
    }

    public double getecan() {
        return ecan;
    }

    private double scgcab = 1;

    // Scaling factor for Max. conductance of Ca background
    public synchronized void setscgcab(double sscgcab) {
        scgcab = sscgcab;
    }

    public double getscgcab() {
        return scgcab;
    }

    /* Na Background Current */
    private double inab; // Na background current (uA/uF)

    public synchronized void setinab(double sinab) {
        inab = sinab;
    }

    public double getinab() {
        return inab;
    }

    private double gnab; // Max. conductance of Na background (mS/uF)

    public double getgnab() {
        return gnab;
    }

    private double enan; // Nernst potential for Na (mV)

    public synchronized void setenan(double senan) {
        enan = senan;
    }

    public double getenan() {
        return enan;
    }

    private double scgnab = 1;

    // Scaling factor for Max. conductance of Na background
    public synchronized void setscgnab(double sscgnab) {
        scgnab = sscgnab;
    }

    public double getscgnab() {
        return scgnab;
    }

    /* ** CONSTRUCTOR ** */
    public LuoRudyModel_II_2000() {
        super();
        getCellGeometry();
    }

    public void setCurrents() {
        currentINa();
        currentICal();
        currentICat();
        currentIKr();
        currentIKs();
        currentIKi();
        currentIKp();
        currentINaCa();
        currentINaK();
        /* This current is only used when modeling pathological conditions */
        //currentINasCa ();
        /* This current is only used when modeling pathological conditions */
        currentIKatp();
        currentIpCa();
        currentICab();
        currentINab();
        currentIt();
    }

    public void setConcentrations() {
        concentrationNai();
        concentrationKi();
        concentrationNsr();
        concentrationJsr();
        concentrationItr();
        concentrationCai();
    }

    public void doSimulate() {
        setCurrents();
        setConcentrations();

        vnew = v - it * dt;
        dvdtnew = (vnew - v) / dt;

        getCaRelease();

        t = t + dt;
        v = vnew;
        dvdt = dvdtnew;
    }

    /** ***************************************************** */
    /**
     * Cell Geometry
     */
    protected void getCellGeometry() {
        vcell = 1000 * Math.PI * a * a * l;
        ageo = 2 * Math.PI * a * a + 2 * Math.PI * a * l;
        acap = ageo * 2;
        vmyo = vcell * 0.68;
        vmito = vcell * 0.26;
        vsr = vcell * 0.06;
        vnsr = vcell * 0.0552;
        vjsr = vcell * 0.0048;
        vcleft = vcell * 0.12 / 0.88;
    }

    /**
     * This is where CICR is calculated, following dv/dtmax cai is measured, 2
     * ms later cai is measured again, the difference of the two determines the
     * amount of SR release
     */
    protected void getCaRelease() {
        if ((v > -35.0) && (dvdtnew < dvdt) && (dvdtmax == false)) {
            caitzero = cai + cmdn + trpn + dcai + csqn + jsr + nsr;
            dvdtmax = true;
            tcount = 0;
        }

        if ((tcount >= 2) && (tcount < (2 + dt))) {
            caiatwo = cai + cmdn + trpn + dcai + csqn + jsr + nsr;
            dcaitwo = caiatwo - caitzero;
            if (dcaitwo >= dcaith) {
                grelbarcicr = gmaxrel * (dcaitwo - 0.00018)
                        / (kmrel + dcaitwo - 0.00018);
                tcicr = 0;
            } else if (dcaitwo < 0) {
                grelbarcicr = 0;
            } else if ((dcaitwo >= 0) && (dcaitwo < dcaith)) {
                dcai = dcaitwo * 1000.0;
                grelbarcicr = gmaxrel
                        * (7.5 * (Math.pow(dcai, 3.0)) - (Math.pow(dcai, 2.0)) + 0.1 * (dcai));
            }
        }

        tcount = tcount + dt;

        /* This is the function that determines if a spontaneous release occurs */

        if ((csqn >= csqnth) && (tjsrol > 50.0)) {
            grelbarjsrol = 4.0;
            tjsrol = 0.0;
        }
    }

    /* Functions that describe the currents begin here */

    /**
     * Calculates Fast Na Current
     */
    protected void currentINa() {
        gna = 16.0 * scgna;
        ena = ((R * temp) / frdy) * Math.log(nao / nai);

        if (v != -47.13) {
            am = 0.32 * (v + 47.13) / (1 - Math.exp(-0.1 * (v + 47.13)));
        } else {
            am = 3.2;
        }
        bm = 0.08 * Math.exp(-v / 11.0);

        if (v < -40.0) {
            ah = 0.135 * Math.exp((80.0 + v) / -6.8);
            bh = 3.56 * Math.exp(0.079 * v) + 310000.0 * Math.exp(0.35 * v);
            aj = (-127140.0 * Math.exp(0.2444 * v) - 0.00003474 * Math
                    .exp(-0.04391 * v))
                    * (v + 37.78) / (1.0 + Math.exp(0.311 * (v + 79.23)));
            bj = (0.1212 * Math.exp(-0.01052 * v))
                    / (1.0 + Math.exp(-0.1378 * (v + 40.14)));
        } else {
            ah = 0.0;
            bh = 1.0 / (0.13 * (1.0 + Math.exp((v + 10.66) / -11.1)));
            aj = 0.0;
            bj = (0.3 * Math.exp(-0.0000002535 * v))
                    / (1.0 + Math.exp(-0.1 * (v + 32.0)));
        }

        if (t == 0) // Steady State Values
        {
            h = ah / (ah + bh);
            j = aj / (aj + bj);
            m = am / (am + bm);
        } else {
            h = ah / (ah + bh) - ((ah / (ah + bh)) - h)
                    * Math.exp(-dt / (1 / (ah + bh)));
            j = aj / (aj + bj) - ((aj / (aj + bj)) - j)
                    * Math.exp(-dt / (1 / (aj + bj)));
            m = am / (am + bm) - ((am / (am + bm)) - m)
                    * Math.exp(-dt / (1 / (am + bm)));
        }

        ina = gna * m * m * m * h * j * (v - ena);
    }

    /**
     * Calculates Currents through L-Type Ca Channel
     */
    protected void currentICal() {
        dss = 1.0 / (1.0 + Math.exp(-(v + 10.0) / 6.24));
        if (v != -10.0) {
            taud = dss * (1.0 - Math.exp(-(v + 10.0) / 6.24))
                    / (0.035 * (v + 10.0));
        } else {
            taud = dss * (1 / 6.24) / 0.035;
        }
        ad = dss / taud;
        bd = (1 - dss) / taud;

        fss = (1.0 / (1.0 + Math.exp((v + 32.0) / 8.0)))
                + (0.6 / (1.0 + Math.exp((50.0 - v) / 20.0)));
        tauf = 1.0 / (0.0197 * Math.exp(-Math.pow(0.0337 * (v + 10.0), 2.0)) + 0.02);
        af = fss / tauf;
        bf = (1 - fss) / tauf;

        if (t == 0) // Steady State Values
        {
            d = ad / (ad + bd);
            f = af / (af + bf);
        } else {
            d = ad / (ad + bd) - ((ad / (ad + bd)) - d)
                    * Math.exp(-dt / (1 / (ad + bd)));
            f = af / (af + bf) - ((af / (af + bf)) - f)
                    * Math.exp(-dt / (1 / (af + bf)));
        }

        if (v != 0) {
            ibarca = pca
                    * zca
                    * zca
                    * ((v * frdy * frdy) / (R * temp))
                    * ((gacai * cai * Math.exp((zca * v * frdy) / (R * temp)) - gacao
                            * cao) / (Math.exp((zca * v * frdy) / (R * temp)) - 1.0));
        } else {
            ibarca = pca * 2 * frdy * (gacai * cai - gacao * cao);
        }
        if (v != 0) {
            ibarna = pna
                    * zna
                    * zna
                    * ((v * frdy * frdy) / (R * temp))
                    * ((ganai * nai * Math.exp((zna * v * frdy) / (R * temp)) - ganao
                            * nao) / (Math.exp((zna * v * frdy) / (R * temp)) - 1.0));
        } else {
            ibarna = pna * frdy * (ganai * nai - ganao * nao);
        }
        if (v != 0) {
            ibark = pk
                    * zk
                    * zk
                    * ((v * frdy * frdy) / (R * temp))
                    * ((gaki * ki * Math.exp((zk * v * frdy) / (R * temp)) - gako
                            * ko) / (Math.exp((zk * v * frdy) / (R * temp)) - 1.0));
        } else {
            ibark = pk * frdy * (gaki * ki - gako * ko);
        }

        fca = 1.0 / (1.0 + cai / kmca);

        ilca = d * f * fca * ibarca * scgcal;
        ilcana = d * f * fca * ibarna * scgcal;
        ilcak = d * f * fca * ibark * scgcal;

        ilcatot = ilca + ilcana + ilcak;
    }

    /**
     * Calculates Currents through T-Type Ca Channel
     */
    protected void currentICat() {
        bss = 1.0 / (1.0 + Math.exp(-(v + 14.0) / 10.8));
        taub = 3.7 + 6.1 / (1.0 + Math.exp((v + 25.0) / 4.5));
        ab = bss / taub;
        bb = (1 - bss) / taub;

        gss = 1.0 / (1.0 + Math.exp((v + 60.0) / 5.6));
        if (v <= 0)
            taug = -0.875 * v + 12.0;
        else
            taug = 12.0;
        ag = gss / taug;
        bg = (1 - gss) / taug;

        if (t == 0) // Steady State Values
        {
            b = ab / (ab + bb);
            g = ag / (ag + bg);
        } else {
            b = ab / (ab + bb) - ((ab / (ab + bb)) - b)
                    * Math.exp(-dt / (1 / (ab + bb)));
            g = ag / (ag + bg) - ((ag / (ag + bg)) - g)
                    * Math.exp(-dt / (1 / (ag + bg)));
        }

        gcat = 0.05 * scgcat;
        eca = (R * temp / (2.0 * frdy)) * Math.log(cao / cai);

        icat = gcat * b * b * g * (v - eca);
    }

    /**
     * Calculates Rapidly Activating K Current
     */
    protected void currentIKr() {
        gkr = (0.02614 * Math.sqrt(ko / 5.4)) * scgkr;
        ekr = ((R * temp) / frdy) * Math.log(ko / ki);

        xrss = 1.0 / (1.0 + Math.exp(-(v + 21.5) / 7.5));
        if (v != -14.2 || v != -38.9) {
            tauxr = 1.0 / (0.00138 * (v + 14.2)
                    / (1.0 - Math.exp(-0.123 * (v + 14.2))) + 0.00061
                    * (v + 38.9) / (Math.exp(0.145 * (v + 38.9)) - 1.0));
        } else {
            if (v == -14.2) {
                tauxr = 1.0 / (0.00138 / 0.123 + 0.00061 * (v + 38.9)
                        / (Math.exp(0.145 * (v + 38.9)) - 1.0));
            } else {
                if (v == -38.9) {
                    tauxr = 1.0 / (0.00138 * (v + 14.2)
                            / (1.0 - Math.exp(-0.123 * (v + 14.2))) + 0.00061 / 0.145);
                }
            }
        }

        axr = xrss / tauxr;
        bxr = (1 - xrss) / tauxr;

        if (t == 0) // Steady State Values
        {
            xr = axr / (axr + bxr);
        } else {
            xr = axr / (axr + bxr) - ((axr / (axr + bxr)) - xr)
                    * Math.exp(-dt / (1 / (axr + bxr)));
        }

        r = 1.0 / (1.0 + Math.exp((v + 9.0) / 22.4));

        ikr = gkr * xr * r * (v - ekr);
    }

    /**
     * Calculates Slowly Activating K Current
     */
    protected void currentIKs() {
        gks = (0.433 * (1.0 + 0.6 / (1.0 + Math.pow(0.000038 / cai, 1.4))))
                * scgks;
        eks = ((R * temp) / frdy)
                * Math.log((ko + prnak * nao) / (ki + prnak * nai));

        xsss = 1.0 / (1.0 + Math.exp(-(v - 1.5) / 16.7));

        if (v != -30.0) {
            tauxs1 = 1.0 / (0.0000719 * (v + 30.0)
                    / (1.0 - Math.exp(-0.148 * (v + 30.0))) + 0.000131
                    * (v + 30.0) / (Math.exp(0.0687 * (v + 30.0)) - 1.0));
        } else {
            tauxs1 = 1.0 / 0.0000719 / 0.148;
        }
        tauxs2 = tauxs1 * 4;

        axs1 = xsss / tauxs1;
        axs2 = xsss / tauxs2;
        bxs1 = (1 - xsss) / tauxs1;
        bxs2 = (1 - xsss) / tauxs2;

        if (t == 0) // Steady State Values
        {
            xs1 = axs1 / (axs1 + bxs1);
            xs2 = axs2 / (axs2 + bxs2);
        } else {
            xs1 = axs1 / (axs1 + bxs1) - ((axs1 / (axs1 + bxs1)) - xs1)
                    * Math.exp(-dt / (1 / (axs1 + bxs1)));
            xs2 = axs2 / (axs2 + bxs2) - ((axs2 / (axs2 + bxs2)) - xs2)
                    * Math.exp(-dt / (1 / (axs2 + bxs2)));
        }

        iks = gks * xs1 * xs2 * (v - eks);
    }

    /**
     * Calculates Time-Independent K Current
     */
    protected void currentIKi() {
        gki = (0.75 * (Math.sqrt(ko / 5.4))) * scgki;
        eki = ((R * temp) / frdy) * Math.log(ko / ki);

        aki = 1.02 / (1.0 + Math.exp(0.2385 * (v - eki - 59.215)));
        bki = (0.49124 * Math.exp(0.08032 * (v - eki + 5.476)) + Math
                .exp(0.06175 * (v - eki - 594.31)))
                / (1.0 + Math.exp(-0.5143 * (v - eki + 4.753)));

        kin = aki / (aki + bki);

        iki = gki * kin * (v - eki);
    }

    /**
     * Calculates Plateau K Current
     */
    protected void currentIKp() {
        gkp = 0.00552 * scgkp;
        ekp = eki;

        kp = 1.0 / (1.0 + Math.exp((7.488 - v) / 5.98));

        ikp = gkp * kp * (v - ekp);
    }

    /**
     * Calculates Na-Ca Exchanger Current
     */
    protected void currentINaCa() {
        inaca = knaca
                * (1.0 / (Math.pow(kmnae, 3.0) + Math.pow(nao, 3.0)))
                * (1.0 / (kmcae + cao))
                * (1.0 / (1.0 + ksat
                        * Math.exp(((nu - 1.0) * v * frdy) / (R * temp))))
                * (Math.pow(nai, 3.0) * cao
                        * Math.exp((nu * v * frdy) / (R * temp)) - Math.pow(
                        nao, 3.0)
                        * cai * Math.exp(((nu - 1.0) * v * frdy) / (R * temp)));
    }

    /**
     * Calculates Na-K Pump Current
     */
    protected void currentINaK() {
        sigma = (Math.exp(nao / 67.3) - 1.0) / 7.0;

        fnak = 1.0 / (1.0 + 0.1245 * Math.exp((-0.1 * v * frdy) / (R * temp)) + 0.0365
                * sigma * Math.exp((-v * frdy) / (R * temp)));

        inak = ibarnak * fnak * (1.0 / (1.0 + Math.pow(kmnai / nai, 1.5)))
                * (ko / (ko + kmko));
    }

    /**
     * Calculates Non-Specific Ca-Activated Current
     */
    protected void currentINasCa() {
        if (v != 0.0) {
            ibarnsna = pnsca
                    * zna
                    * zna
                    * ((v * frdy * frdy) / (R * temp))
                    * ((ganai * nai * Math.exp((zna * v * frdy) / (R * temp)) - ganao
                            * nao) / (Math.exp((zna * v * frdy) / (R * temp)) - 1.0));
        } else {
            ibarnsna = pnsca * frdy * (ganai * nai - ganao * nao);
        }
        if (v != 0.0) {
            ibarnsk = pnsca
                    * zk
                    * zk
                    * ((v * frdy * frdy) / (R * temp))
                    * ((gaki * ki * Math.exp((zk * v * frdy) / (R * temp)) - gako
                            * ko) / (Math.exp((zk * v * frdy) / (R * temp)) - 1.0));
        } else {
            ibarnsk = pnsca * frdy * (gaki * ki - gako * ko);
        }

        insna = ibarnsna / (1.0 + Math.pow(kmnsca / cai, 3.0));
        insk = ibarnsk / (1.0 + Math.pow(kmnsca / cai, 3.0));
    }

    /**
     * Calculates ATP-sensitive Potassium Current
     * <p>
     * Ferrero et al., Circ. Res. 1997;79:208-21.
     */
    protected void currentIKatp() {
        //channel density
        double satp = 3.8;
        //maximum channel probability
        double p0 = 0.91;

        //unitary conductance in the absence of ATP
        //double gamma0 = 35.375 * (Math.pow((ko / 5.4), 0.24));
        //modified to account for the steeper Ko-dependence
        double gamma0 = 35.375 * (Math.pow((ko / 5.4), 0.54));

        //rectification caused by Mg
        double khmg = (0.65 / Math.sqrt(ko + 5.0))
                * Math.exp(-(0.64 * frdy * v) / (R * temp));
        double fm = 1.0 / (1.0 + (mgi / khmg));

        //rectification caused by Na
        double khna = 25.9 * Math.exp(-(0.35 * frdy * v) / (R * temp));
        double fn = 1.0 / (1.0 + Math.pow((nai / khna), 2.0));

        //Temperature effect
        double ft = 1.0; //Math.pow(1.3, ((temp - 36.0) / 10.0));

        //unitary conductance
        double g0 = gamma0 * fm * fn * ft;

        //fraction of open KATP channels
        double kmatp = 35.8 + 17.9 * Math.pow(adpi, 0.256);
        double hatp = 1.3 + 0.74 * Math.exp(-0.09 * adpi);

        fatp = 1.0 / (1.0 + Math.pow((atpi / kmatp), hatp));

        gbarkatp = satp * scgbarkatp;

        //From Shawn, Rudy, Cardiovasc. Res. 1997;35:256-72
        //patp = 1 / (1 + (Math.pow(3 / .25, 2)));
        //ikatp = gbarkatp * patp * (Math.pow((ko / 4), 0.24)) * (v - eki);

        ikatp = gbarkatp * g0 * p0 * fatp * (v - eki);
    }

    /**
     * Calculates Sarcolemmal Ca Pump
     */
    protected void currentIpCa() {
        ipca = (ibarpca * cai) / (kmpca + cai);
    }

    /**
     * Calculates Ca Background Current
     */
    protected void currentICab() {
        gcab = 0.003016 * scgcab;
        ecan = ((R * temp) / (2.0 * frdy)) * Math.log(cao / cai);

        icab = gcab * (v - ecan);
    }

    /**
     * Calculates Na Background Current
     */
    protected void currentINab() {
        gnab = 0.00141 * scgnab;
        enan = ena;

        inab = gnab * (v - enan);
    }

    /**
     * Calculates Total Na Current
     */
    protected double getTotalINa() {
        return getina() + getinab() + getilcana() + getinsna() + 3 * getinak()
                + 3 * getinaca();
    }

    /**
     * Calculates Total K Current
     */
    protected double getTotalIK() {
        return getikr() + getiks() + getiki() + getikp() + getilcak()
                + getinsk() + getikatp() - 2 * getinak();
    }

    /**
     * Calculates Total Ca Current
     */
    protected double getTotalICa() {
        return getilca() + geticab() + getipca() - 2 * getinaca() + geticat();
    }

    /**
     * Calculates Total Current
     *
     * Total sum of currents is calculated here, if the time is between stimtime =
     * 0 and stimtime = tstim, a stimulus is applied
     */
    public void currentIt() {
        naiont = getTotalINa();
        kiont = getTotalIK();
        caiont = getTotalICa();

        if ((t >= tstim) && (t < (tstim + dt))) {
            stimtime = 0;
            tstim = tstim + bcl;
            dvdtmax = false;
        }

        if ((stimtime >= 0) && (stimtime < stl)) {
            it = st + naiont + kiont + caiont;
        } else {
            it = naiont + kiont + caiont;
        }

        stimtime = stimtime + dt;
    }

    /* Functions that calculate intracellular ion concentrations */

    /**
     * Calculates new myoplasmic Na ion concentration
     */
    protected void concentrationNai() {
        dnai = -dt * (naiont * acap) / (vmyo * zna * frdy);
        nai = dnai + nai;
    }

    /**
     * Calculates new myoplasmic K ion concentration
     */
    protected void concentrationKi() {
        dki = -dt * (kiont * acap) / (vmyo * zk * frdy);
        ki = dki + ki;
    }

    /**
     * Calculates new NSR Ca ion concentration
     */
    protected void concentrationNsr() {
        kleak = iupbar / nsrbar;
        ileak = kleak * nsr;
        iup = iupbar * cai / (cai + kmup);

        dnsr = dt * (iup - ileak - itr * vjsr / vnsr);
        nsr = nsr + dnsr;
    }

    /**
     * Calculates new JSR Ca ion concentration
     */
    protected void concentrationJsr() {
        grelcicr = grelbarcicr * (1.0 - Math.exp(-tcicr / tauon))
                * Math.exp(-tcicr / tauoff);
        irelcicr = grelcicr * (jsr - cai);
        tcicr = tcicr + dt;

        greljsrol = grelbarjsrol * (1.0 - Math.exp(-tjsrol / tauon))
                * Math.exp(-tjsrol / tauoff);
        ireljsrol = greljsrol * (jsr - cai);
        tjsrol = tjsrol + dt;

        csqn = csqnbar * (jsr / (jsr + kmcsqn));
        djsr = dt * (itr - irelcicr - ireljsrol);
        bjsr = csqnbar - csqn - djsr - jsr + kmcsqn;
        cjsr = kmcsqn * (csqn + djsr + jsr);

        jsr = (Math.sqrt(bjsr * bjsr + 4.0 * cjsr) - bjsr) / 2.0;
    }

    /**
     * Calculates Translocation of Ca from NSR to JSR
     */
    protected void concentrationItr() {
        itr = (nsr - jsr) / tautr;
    }

    /**
     * Calculates new myoplasmic Ca ion concentration
     */
    protected void concentrationCai() {
        dcai = -dt
                * (((caiont * acap) / (vmyo * zca * frdy))
                        + ((iup - ileak) * vnsr / vmyo)
                        - (irelcicr * vjsr / vmyo) - (ireljsrol * vjsr / vmyo));
        trpn = trpnbar * (cai / (cai + kmtrpn));
        cmdn = cmdnbar * (cai / (cai + kmcmdn));

        catotal = trpn + cmdn + dcai + cai;
        bmyo = cmdnbar + trpnbar - catotal + kmtrpn + kmcmdn;
        cmyo = (kmcmdn * kmtrpn) - (catotal * (kmtrpn + kmcmdn))
                + (trpnbar * kmcmdn) + (cmdnbar * kmtrpn);
        dmyo = -kmtrpn * kmcmdn * catotal;

        gpig = Math.sqrt(bmyo * bmyo - 3.0 * cmyo);

        cai = (2.0 * gpig / 3.0)
                * Math
                        .cos(Math.acos((9.0 * bmyo * cmyo - 2.0 * bmyo * bmyo
                                * bmyo - 27.0 * dmyo)
                                / (2.0 * Math.pow((bmyo * bmyo - 3.0 * cmyo),
                                        1.5))) / 3.0) - (bmyo / 3.0);
    }
}